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We consider dimensional reduction techniques for the Liouville-von Neumann equation for the evaluation of 
the expectation values in a mixed quantum system. In applications such as nuclear spin dynamics the main 
goal for simulations is being able to simulate a system with as many spins as possible, for this reason it is 
very important to have an efficient method that scales well with respect to particle numbers. We describe 
several existing methods that have appeared in the literature, pointing out their limitations particularly in 
the setting of large systems. We introduce a method for direct computation of expectations via Chebyshev 
polynomials (DEC) based on evaluation of a trace formula combined with expansion in modified Chebyshev 
polynomials. This reduction is highly efficient and does not destroy any information. We demonstrate the 
practical application of the scheme for a nuclear spin system and compare with several alternatives, focusing 
on the performance of the various methods with increasing system dimension. Our method may be applied 
to autonomous quantum problems where the desired outcome of quantum simulation, rather than being a 
full description of the system dynamics, is only the expectation value of some given observable. 
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I. INTRODUCTION 

The motion of a quantum system of particles is de- 
scribed by the Schrodinger equation : 



gjjg) 



= -iH|*>, 



(1) 



where H is the Hamiltonian and the wave function of 
the system, and we have chosen physical quantities such 
that h — 1. When H is time independent it is possible 
to have an exact representation of the solution: 
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(2) 



From a numerical point of view finding a solution of Q is 
a challenging task, due to the cost of computing the op- 
erator exponential e~ lHt . Many numerical methods have 
been proposed over the years to solve ([2]), from polyno- 
mial expansion of the exponential*—, to projection on 
Krylov subspaces^. When the wave function \^f) de- 
pends on the particle positions q and momenta p it is 
very common to adopt a splitting method for the Hamil- 
tonian, H — T + V being V diagonal and T the kinetic 
term which is diagonal in Fourier space. With this ap- 
proach the main issue is the evaluation of a fast Fourier 
transform (FFT) to switch between the x-grid and the p- 
grid; the cost of this transform is of order 0(N log N)^. 
However when the wave function depends on other vari- 
ables, like the spin of the particles, it is not straightfor- 
ward to apply a splitting method. For a general survey 
of these approaches see Refs. and|||. 

It is important to remark that all of these methods 
suffer when the dimension of the space needed to repre- 
sent the wave function $ accurately becomes large. In 
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many practical cases, as for example in spin dynamics, 
the complexity of the computational problem grows ex- 
ponentially with the number of atoms. For pure spin 
dynamics "J may be represented by a vector in a space of 
dimension (2/+ 1)™, where n is the number of particles 
with spin/ present. 

This paper is organized as follows: in §2 we introduce 
the Liouville-von Neumann equation that describes the 
dynamics of the density matrix, we also clarify the issues 
that complicate its numerical solution. In §3 we discuss 
two of the main techniques that have been applied in 
quantum simulation for the evaluation of the exponential 
that appears in the solution of the Schrodinger equation, 
specifically Krylov expansion via Lanczos-Arnoldi itera- 
tions and the expansion in Chebyshev polynomials. We 
also present a more recent method based on a different 
evaluation of the Krylov subspace expansion called Zero 
Track Elimination (ZTE) 23 . This method has been de- 
veloped for simulations of nuclear spin dynamics within 
the Nuclear Magnetic Resonance (NMR) community. We 
give what we believe to be the first mathematical discus- 
sion of the method, exploring its limitations and clarify- 
ing situations in which it works. 

In §4 we introduce a new method we have developed: 
the Direct Expectation values via Chebyshev (DEC). 
This is the main result of this paper. In this section we 
also discuss the technical details of the implementation 
of the DEC method. 

In §5 we provide a performance comparison of DEC 
with three separate alternatives: the classic Krylov ex- 
pansion based on Lanczos algorithm, with the new ZTE 
method and with the Chebyshev expansion. The sample 
system considered is a pairwise interacting nuclear spin 
system; even if it is a simple system it contains the key 
features of the systems of interest to the NMR commu- 
nity. 
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II. THE LIOUVILLE VON-NEUMANN EQUATION 

When the sample of interest is composed of a large 
number N of identical systems each described by a wave 
function \^), j = 1, . . . , N, we may introduce a statisti- 
cal operator called density operator g defined by: 

N 

q = J2\V)(V\, (3) 

J=l 

which evolves according to the Liouville-von Neumann 
equation: 

^=-i[H,g}. (4) 

If g is expanded using a (finite) approximate basis set 
{\(fi) ■ ■ ■ \<p n )}> © ma y be viewed as an ordinary differ- 
ential equation in a matrix argument. In the time inde- 
pendent case the solution for ^ may be written: 

g(t) = e- lHt g Q e im . (5) 

It is possible and sometimes preferable to rewrite ([5]) 
by use of the Liouvillian L = Id <g> H — H ® Id, where Id is 
the identity matrix, allowing us to recast g as a vector: 

g(t) = e- lLt g n . (6) 

The main issue is then to evaluate the exponential of the 
matrix L. For a large system and using a typical basis 
set, H (and consequently L) will be a structured sparse 
matrix. 

While the evolution of the system is described by the 
density matrix, the outputs we are typically interested 
in obtaining from quantum simulations are the expecta- 
tions of observables, these being the only quantities we 
can compare with the experiments. In the density ma- 
trix formalism the expectation value of an observable Q, 
associated with an operator Q is written as: 

(Q(t)) = Trace{ e (t)Q}. (7) 

Whereas in general quantum simulations the equation of 
motion is solved for a quantity, g, which has dimension 
n x n, the types of outputs we are actually interested in 
are typically one dimensional objects ([7]). In this paper 
we exploit this fact and design an algorithm that com- 
putes (almost) directly the evolution of the expectation 
value (Jjj , instead of the evolution of the density matrix 
([5]). This approach does not lose any information of the 
original system (the only errors arise due to truncation), 
but at the same time the method provides a powerful 
computational tool, with potential dramatic reduction 
in the computational cost, especially when dealing with 
large matrices. The main idea of this approach is to ex- 
ploit features of a Chcbyshev expansion for the matrix 
exponential in (J6j) . 

Several methods to evaluate ^ are discussed in a re- 
cent monograph-; however our proposed Direct evalua- 
tion of the Expectation values via Chebyshev polynomials 



(DEC) method is different because it does not solve the 
evolution for the density matrix. Instead it exploits the 
trace evaluation in ([7]) and with the evaluation of just 
one Chebyshev expansion it allows the solution of (J7J at 
any time. DEC can be extremely powerful when we are 
only interested in the expectation values; if instead it is 
necessary to evaluate the evolution of the density matrix 
itself ©, then a traditional approach might be the best 
choice. 



III. EXISTING METHODS 

A. Krylov expansion via Lanczos-Arnoldi 

In the past decades many methods have been 
developed for numerical evaluation of the matrix 
exponential, 9 . However for large sparse matrices, the 
methods usually applied are those based on expansion 
in Krylov subspaceiSr— . The main idea is to project © 
onto the subspace: 

K m (L, g) = span{p , Lg , L 2 g , . . . , L m g }. (8) 

To get a suitable basis set of ([5]) we may use the Lanczos 
algorithm 1 ^, as L is Hermitian. The Lanczos method 
is an iterative method, a very desirable feature in the 
context of large sparse matrices. For specific details see 
Ref^i 6 .. 

An approximation for g(t) is: 

g(t) = e- lLt g Q a \\g \\V m e- iTmt ex, (9) 

where T m and V m come from the Lanczos algorithm. 
The Lanczos algorithm provides an orthonormal basis set 
V m for the Krylov subspace K m (L, go) via a three-term 
recursio n 14 ' 16 : 

Pj+iqj+l = Lqj - a-jq-j - foqj-i qi = g . (10) 

T m is a tridiagonal matrix of size mxm, and e\ is the first 
vector of the canonical basis of size n. This technique is 
very powerful for short time simulations, because with 
few iterations m it is possible to have remarkably good 
approximations, but for longer times larger Krylov sub- 
spaces would be needed to stay close to the real solution. 
On the other hand if we do not consider enough terms in 
the Lanczos algorithm for longer times, (J9j) is no longer 
a reliable approximation. 

It is possible to set a stopping criterion for the Lanczos 
iterations^ 1 -; for a given t we can find m such that: 

ipm]m+l,m||e %tTm || m ,l < £. (11) 

In our experiments the best way to implement a Krylov 
expansion was to evaluate each step: 

Qn+i = e- lLdt g n « WQnWV^e-^ex. (12) 

In this way with less than 10 iterations of Lanczos per 
step it is possible to have a fast and accurate benchmark. 
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The obvious drawback is that no information passes from 
step n to step n + 1. However in our numerical tests 
the use of a longer timestep that would allow a common 
Krylov subset K m (L,g) for more than one step g n was 
not preferable as more iterations of (ITOl) were then re- 
quired. As the equation (fTT|) involves the evaluation of 
the exponential of a tridiagonal matrix, when m becomes 
large, this operation may become a serious bottleneck for 
the whole simulation. 



B. The Chebyshev expansion 

Another method that has been widely applied for solv- 
ing ([!} is the expansion of the exponential in Chebyshev 
polynomials^. 

The preliminary step of this method is to rescale the 
matrix within the interval [—1, 1], as outside this inter- 
val the Chebyshev polynomials grow rapidly, and the ex- 
pansion becomes unstable; to do that we need to eval- 
uate the two extremes of the spectrum of L. In or- 
der to obtain extreme values we propose, as already 
mentioned in the literature^, to perform a few steps 
of Lanczos iteration, as this provides a good approxi- 
mation for the extreme eigenvalues, for small compu- 
tational cost. If we define these two values as a and 
/3, and assume /3 < <r(L) < a, we may rewrite L as 
L = (SId-L s D), where D = (a-/3)/2, S = (a + /3)/2, 
and —1 < u(L s ) < 1. We may then expand the exponen- 
tial of L s in the Chebyshev polynomials and we arrive at 
the following equation for g: 

e {t) = e - iLt g a w e- its ( c k (t D )T k (L s )g \ , (13) 
V fc=o / 

with to = Dt. Both c k (to) and T k (L s ) can be calculated 
iteratively: 

c fc (t) = {2-5 kfi )(-i) k J k (t), (14) 



start the backward iteration process from m s t a rt = n + r, 
where n is the actual order of the function we are inter- 
ested on and r is some small expansion. In this case we 
need to know already from an a priori error analysis how 
many iterations need to be performed to to get below the 
threshold e. 

It is possible to prove that for the rescaled Hcrmitian 
matrix L s , when applied to a vector of unit Euclidian 
norm, we have^: 

\\P m -x{tL s )Q Q -e- itL 'Q \\<A{e l -^ 2m ^^-) m , (17) 

2m 

for m > t, where P m (t) is the order m expansion in 
Chebyshev polynomials. This equation indicates that 
there is a superlinear decay of the error when m > t. 
The formula may be derived by examining the asymp- 
totic behaviour of the Bessel functions (|37|) . We may 
then use the relation 4(exp{l — (t /2m) 2 }^) m < e to 
approximate m. 

From a practical point of view the usual way of apply- 
ing Chebyshev is to evaluate^: 

Qn+i = e- lLdt g n ~ P nn ^(dtL)g n , (18) 

where: 

™max 

Pn n , ax = e~ ldts ]T c k (dtD)T k (L s ). (19) 

fc=0 

m can be evaluated either from (IT7|) or directly checking 
the convergence (f3"5)) . 

To avoid numerical instabilities coming from the iter- 
ative formula for the Bessel functions it is also possible 
to get P m for a given m via the Clenshaw Algorithm^^: 

d k = c k g n +2L s d k+ i-d k+2 , k = m—l,m—2, . . . , 0, (20) 

with initial values d TO +i = d m = 0, and P rn (dtL)go = 
d - d 2 . 



T k+1 (x) = 2T k {x)x - T fc _i(a:), (15) 

with initial values Tq(x) = Id, Ti(x) = x. J k (t) is the 
fc-th Bessel function of the first kind. 

The Bessel Functions of the First Kind of integer order 
may be evaluated directly by using a three-term recur- 
rence relation 

2n 

J n +i(t) = —J n - J n -l- (!6) 

It is well known that (|T5|) becomes numerically unstable 
for n > t, 1 -. To improve the method, we may exploit 
the linear nature of the iterative algorithm. It is pos- 
sible to use Miller's algorithm, and to solve an inverted 
form of (fTBT) . i.e. to compute J„_i given J n , J n+1 21. 
When using Miller's Algorithm it is suggested to expand 
the number of terms (providing a sort of buffer), i.e. to 



C. The Zero-Track-Elimination method 

Nuclear Magnetic Resonance (NMR) is a spectroscopy 
technique that exploits the interaction between nuclear 
spins and electromagnetic fields in order to analyse the 
samples. The temporal evolution of such a system is de- 
scribed via a density matrix that has size (2,1 + 1/2)" 
where / is the spin and n the number of nuclei. The 
exponential growth of the size of g with respect to n im- 
pedes the use of simulations when dealing with systems 
involving more than few (5-6) spins. Many attempts have 
been made to solve this (see Refs. [13 and [23| for recent 
approaches), even using Chebyshev polynomials^. These 
algorithms have been developed to simulate both liquid 
systems, where the Hamiltonian is generally time inde- 
pendent, and for Solid-State NMR. In the last case the 
Hamiltonian is time dependent due to the non averaging 
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out of anisotropic interactions during the motion of the 
sample, and as previously remarked DEC is not applica- 
ble in this case. 

Recently also a new method for model reduction in the 
simulation of large spin systems, so called Zero Track 
Elimination (ZTE), has been presented^. This tech- 
nique is based on the idea of pruning out the elements 
of g(t) which do not belong to K (L, qq). In practice, the 
initial density matrix go is a sparse vector. 

In order to reduce the steps needed to evolve the full 
system, we monitor the elements of g(t) that stay be- 
low a chosen threshold £ during this first evolution steps 
and introduce structural zeros based on these observa- 
tions. The evolution is then performed in this reduced 
state space (gz,Lz)- The idea is extremely appealing, 
as once the propagator for L z is evaluated all the sub- 
sequent steps have the cost of a reduced matrix- vector, 
and it is possible to use (|4"P|) for the reduced system. 

The initial time length is set as the inverse of the 
largest Larmor frequency. The Larmor frequency is a 
given quantity for each element that depends on the phys- 
ical property of the nucleus and is the frequency of res- 
onance for a non-interacting spin: oj® — ~jjB, where 
— 7 3 is the gyromagnetic ration of the nucleus and B the 
applied magnetic field. 

The time-length of the initial check is then: 



ti ar minj{\uj®\} 

where w° is the Larmor frequency of the j-th spin. The 
following theorem given in^ assures that it is possible 
to prune out from the evolution those states that remain 
exactly during the first few time steps. It is possible to 
prove that for a state \l) there holds the equation: 

(l\e- iLt \g Q ) =0,t€ [0, St] (l\er lLt \g Q ) = 0,t£ [0, oo) 

(22) 

To prove (f2"2")) it is sufficient to expand e lLt into a 
Taylor expansion. In fact: 

(l\e-* Lt \go) = (l\ £ L k( ~^\g ) = £ tB-(l\L k Qo ), 

k=0 ' k=Q 

(23) 

can be true for any t G [0,5t] only if (l\L k \g ) = for all 
the k, but this means that (l\e~ lLt \go) = 0,tt G [0, oo). 

From a practical point of view few states (I | obey (|221) , 
but a much higher number of states will stay close to 
during the first j time steps, where j — St/dt. So the 
states (l\ that are pruned out are those for which: 

(l\e- iLt \Q ) < e, t&[0,Si] (24) 

It is claimed^ that the error of such an approximation 
is similar to what would be obtained by considering in 
the Krylov expansion the contributions coming from high 
values of n in L n g . 

The main problem of this approach is the choice of St, 
i.e. the duration of the initial propagation. In fact if we 



look at a full diagonalization of L we see that: 

g(t) = Xe- wt X- l Q , (25) 

so each element of g(t) can be written as a sum of oscil- 
lators vibrating at the different Xj € u(L). 

n N 

g k {t) = ]T X [Km e-^\ H = ]T X bA gl (26) 
j=i i=i 

From (|26p it is clear that to ensure that we are not prun- 
ing out the low frequencies modes we would need at least 
St oc 1/ minj { | A j | } , and this could be different from the 
lowest Larmor frequency min^ { | w° | } as this last is a quan- 
tity related to the non-interacting spins. 

It is possible to have an estimate of the lowest eigen- 
value of L using a technique like the one we used for the 
Chebyshev expansion QUI Bp . The simple choice of St de- 
pending on Xj rather then w° is not enough to ensure the 
validity of (12"4"1). 

Within this framework we can restate (|24[) in a simpler 
form for a one dimension function and prove that: 

Theorem III.l Given a function f(t) = X^j=i e_ J^j 
s.t. \f(t)\ < e for t £ [0, 27r/ minjjAj}], this is not suffi- 
cient to ensure that \ f{t)\ < e for t £ [0, oo). 

Proof To prove (|III.1[) we can check that if the Xj are not 
well separated then it is possible to have a combination of 
similar frequencies that build up on a total frequency that 
is the lowest common multiplier of the initial frequencies. 
In the general form if for one or more g k we have the 
resonance condition: 

g k (t) = ae-^ xt - £ po e ^{X+P)t^ (2?) 
j 

and fij , a > rationals such that a — • f3j . 3 small 
enough Sj s.t. we have that |£>fc(£)l < e, t £ [0,St] but at 
the periodic maximum t > St s.t. \f(t)\ = a+Y] ■ f3j 3> e. 

In the appendix we present one example of such a func- 
tion. 

In order to ensure the validity of (j2~4")l even with St that 
depends on a(L) rather than the Larmor frequencies, we 
need to check the separation of the eigenvalues, and ob- 
viously an analysis of that kind would be as expensive as 
a whole simulation. 

In practice, it might be argued that we are unlikely 
to have a dynamics like (j2"T)l . but especially with systems 
with many spins, a wide variety of behaviors are certainly 
possible in the general situation. 

Focusing on NMR simulations however the reason why 
ZTE performs so well22 lies in the fact that the numerical 
comparison with other methods is not performed on g(t) 
but on the observables ([7]) . In particular for experimental 
reasons the only quantity that can be compared with 
experiments is the free induction decay (FID) signal: 

f{t) = Trace {g(t)I p } . (28) 
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g(t) can be written as combination of Pauli matrices, and 
Ip is the shift up operator: I p = I x + il y - The Fourier 
transform of f(t) gives then the spectrum of the sample, 
where each resonance frequency is revealed by a peak. 
Due to relaxation effects in the experiments the shape 
of the resonance peaks is not a delta function (as would 
result from (|26[) but is a smooth function. To model 
this smoothness also for the simulated data it is common 
practice to evaluate the Fourier transform not of (|28|) but 
for an exponentially decaying function f(t): 

f(t) = e-«*/(t), (29) 

where the parameter £ comes from other fitted data 21 . 
The main effect of (l29l) is exactly to smooth out all 
the low frequency modes that will differentiate g(t) from 
Qz{t). 

There are however some drawbacks: 

• for this method there is no available convergence 
theory; 

• the performance depends strongly from the initial 
condition go, and on H. As expected, in our tests 
the size of the reduced system could change by a 
factor 2 depending on the number of interacting 
spins. The reason for this effect comes from the fact 
that the less sparse L is, the more non-zero states 
will appear within the first steps, thus reducing the 
effectiveness of the reduction L z . 

• Another reason of the strong dependence of ZTE 
with respect to the initial conditions comes from 
(1211) : depending on the Larmor frequencies, and on 
the timestep size the number of evolution steps at 
the beginning can become large, and this being the 
most expensive part of the simulation the influence 
of it on the total computation costs can become 
significant. 

IV. THE DIRECT EXPECTATION VALUES VIA 
CHEBYSHEV 

The common point of all the methods presented in the 
previous section is that they involve the propagation of 
the matrix g, and for this reason they suffer from re- 
quiring that matrix operations (or matrix-vector oper- 
ations) be performed at each step of calculation. Let 
us recall that the only quantities that can be compared 
between quantum simulations and experiments are the 
observables, for an operator Q we have: 

(Q(t)) = Trace{g(t)Q}. (30) 

The first step of DEC is to perform a Chebyshev expan- 
sion as seen in section UlI Bl to get: 

e (t) = e-* Lt g a e-* ts c k (t D )T k (L s )g Q J , (31) 

V fe=o / 



If we insert (l3Tj) into (l3Cfl) we find: 

(Q(t)) = Trace j |V ltS "f] c k {t D )T k {L s )gA Q j . 

'7" . < 32 ) 

By exploiting the linearity of the trace operation we can 

pull out of the trace all time dependent parts, and eval- 
uate once for all the coefficients T k (L s ). In fact we may 
rewrite ([3"2l as: 

™max 

(Q(t)) = e- lts ]T Cfe^Trace {r fc (L s )<)} . (33) 

k=0 

This is the key equation of the DEC method as it 
is possible to store an array of scalar values T k — 
Trace{(Tfc£>o)Q}- All the time dependent terms are just 
scalar values that have to be multiplied by T k to get the 
evolution of Q at any time: 

n maK 

(Q(t)) =e~ us ]T c k (tD)f k . (34) 

k=0 

If more than one observable is required it is still pos- 
sible to use DEC. The only difference with the single 
expectation case is that we need to store different sets of 
Tfc, one for each operator Q. 

A. Stopping Criterion 

The number of terms for the polynomial expansion 
in (Til?)) depends on a prescribed tolerance e, and on 
the time tr>. In other methods based on Chebyshev 
approximation^, the following has been suggested as a 
stopping criterion: 

»1max S.t. \\Cn mw {t D )\\ < E. (35) 

Due to the zeros of the Bessel function J, at fixed time to, 
(|35p may hold for some n, even though the expansion has 
not yet reached the convergence regime; it may happen 
that for rii > n we have that c ni {tu) > c„(i£>). To avoid 
this effect it is enough to use as a stopping criterion a 
combination of two Bessel functions; the cost of such a 
stopping criterion is that at most we need to perform an 
extra iteration step (|14[) . In our numerical tests we have 
used the following: 

(^)|| 2 + ||c„ max (^) 2 || <e. (36) 

The total time r plays a role here, since the larger r the 
more terms (T k ,c k ) will be needed to get \c k \ below the 
threshold e. 



B. Computation of the Expansion 

In order to optimise the number of terms we evaluate, 
but without having to check at each step whether we 
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FIG. 1. Example of few integer order Bessel Functions of the 
first kind. 

have already evaluated enough terms T k , we propose to 
evaluate first (Q(t)), at the final time r, and to store the 
Njobx values of T k . We can prove that: 

Theorem IV. 1 Given a Chebyshev expansion for an ex- 
ponential e~ lLt Qo if i36\) holds for a given time t and 
small enough e then 136\) holds for any time t < r. 

Proof From equation (IT4l) it is clear that c k depends 
on the Bessel functions. If we look at the asymptotic 
behaviour of the Bessel function of fist kind, for any k € 
N, we have that, for k fixed 18 : 

J *W~rWi)G)*' lim ^ ' (37) 

where T{t) is the Euler-r and for n £ Z we have that 
r(n) = (n — 1)!. Equation (|37|l shows that for any k ^ 0, 
in a neighbourhood of t = 0, Jfc(t) is increasing mono- 
tonically with respect to t. This behaviour is maintained 
for the whole interval [0,j' k ] where j' k is the first zero of 
the derivative of J k . It is possible to show (see Ref. [l8|, 
Eq.9.5.2), that k < j' k ; consequently we can say that if 
(fM)) holds for a given n max at r and r < n max , then we 
are in the monotonically increasing region for J nmax and 
J„ max+ i. In this case, equation ([36]) holds also for any 

t <T. 



C. Efficient Implementation 

The cost of DEC is all in the first step. Note that the 
cost of the evaluation of any T k itself is roughly equiva- 
lent to that of a matrix-matrix multiplication, as per the 
iteration T k+1 {L S ) = 2T k (L s ) - T fc _i {HJ. But what is 
actually needed in all our calculations is T^qq. Because 
of the linearity of the iterative expression, we may multi- 
ply To and T\ by go and then use (|T5j) directly on T k go- 
The iterated operation in then just a matrix-vector mul- 
tiplication. 



After all the T k needed have been stored, it is possible 
to get (Q(t)) at any time t £ [0, r] by evaluating: 

(Q(t)) = e- tDt ° J2 c k(ts)f k (38) 
fe=l 

where the c k are evaluated iteratively via (|14p , and m max 
satisfies (l36l) for tD. 

For the details of the algorithm we refer to the Appendix. 
As a final remark we note that if the Hamiltonian is time 
dependent it is not possible to apply DEC, as it is not 
possible anymore to isolate the time dependent part out 
of the trace. 



V. NUMERICAL EXPERIMENTS 

As in many other physical systems, nuclear spin dy- 
namics provides a perfect example to test DEC, because 
the final outcome of the simulations is an observable, the 
free induction decay (FID) signal, and this result is the 
sole important quantity, as it is the only data available 
from experiment. 

As Hamiltonian we assumed a sum of isotropic chem- 
ical shift and the isotropic term of a pair interaction 
called Homonuclear J-couplings, that depends on the in- 
ner product Ij ■ 

n n 

H = -zZ UJ P^+zZ J ^-Ii- (39) 

For the initial density matrix we set go — —ly, that is the 
result of the application of a so called cc-pulse to a sample 
already under the effect of a strong constant magnetic 
field along the z direction™ This is the usual initial 
condition when the acquisition of the signal starts. 

An illustration of the structure of the Liouvillian ma- 
trix is presented in Figj2] The sparsity depends on the 
number of interactions among the spins. In most cases 
the J-coupling interaction matrix J is relatively sparse. 
In our numerical test a strong coupling system, where J 
is dense and each spin interacts with every other spins, 
was simulated. 

Due to the fact that our implementation involves only 
matrix-vector multiplication, techniques developed both 
for structured and unstructured sparse matrices may be 
exploited. 

For comparison of computational costs we tested this 
method with an increasing number of spin particles us- 
ing different methods to evaluate the exponential. In 
particular to examine the error we compared DEC with 
the expm function of MATLAB, that uses a scaling and 
squaring algorithm with Pade' approximation. In this 
way we evaluate once for all U — e~ lLdt where dt is the 
step size of the simulation, and then at each time-step 
we propagate g: 

g n+1 = Ug n . (40) 





5000 10000 15000 

nz = 180096 



FIG. 2. Left: Sparse structure of the Liouvillian (n = 1024, 
nz — 6112) for a system of 5 spins. Right: Structure of the 
Liouvillian (n = 16384, nz = 180096) for a system of 7 spins. 
Both the systems are in a weak coupling condition, each spin 
interacts with approximately half the other spins. 
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FIG. 3. Logarithmic comparison of computational costs for 
Lanczos, Zero Track Pruning (ZTE) and Direct Expectations 
via Chebyshev (DEC), with dt = 0.1 N = 1000. 



It is well known that in terms of computational costs 
this simplistic approach performs poorly, so we compared 
DEC also with the other methods presented in this paper: 
Krylov expansion via Lanczos^, the Chebyshev expan- 
sion, and ZTE22. 

For the Lanczos method we have used the function expv 
of the package EXPOKIT 24 written in textslMATLAB. 
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FIG. 4. Logarithmic comparison of computational costs for 
Lanczos and DEC when simulating for the same number of 
total steps N but with different stepzise dt. N — 1000 in both 
the cases. 



Matrix size 


expm 


expv 


Chebyshev 


ZTE 


Reduced^ 


DEC 


16 


0.08 


0.99 


0.54 


1.06 


8 


0.39 


64 


0.1 


1.41 


0.68 


1.79 


30 


0.60 


256 


0.17 


1.88 


0.95 


2.26 


112 


0.76 


1024 


1.5 


7.87 


2.40 


5.80 


420 


1.14 


4096 


40.84 


29.78 


9.60 


24.51 


1584 


1.93 


16384 




165.06 


102.92 


112.85 


6006 


16.11 


65536 






677.22 






129.93 


262144 












542.98 



a Size of q for the reduced system 

TABLE I. Comparison of computational costs, CPU time in 
seconds, for dt = 0.1, N = 1000. 



these numerical tests we set the tolerance to be e — le—7. 

DEC performs at its best for short time simulations 
(i.e. when total time r is small), so that we do not need 
to evaluate a large number of T& , and when at the same 
time the use of small time step dt is required, as the cost 
for any step after the first is negligible. For instance, 
while for all the other methods the cost of a 1000 step 
simulation with dt — 0.1, is approximately half the cost 
of a simulation of 1000 steps with dt = 0.1, for DEC there 
is a gain of almost an order of magnitude, see Figj4] 



VI. CONCLUSION 



A. Summary of Results 

The error-to-cost (measured in CPU time) diagrams 
are shown in Figure [3] for all the methods described. All 
the numerical tests have been performed on a Dell Pow- 
erEdge 1950 wiht 4GB RAM and a DualCore Intel pro- 
cessors running in 32bit mode. The language used is 
MATLAB. 

It is clear that DEC is almost an order of magnitude 
more efficient than the alternatives. To avoid instabilities 
coming form the evaluation of the Bessel functions in 



In this article we have presented a new method for 
simulation of an observable in a mixed quantum sys- 
tem. By expanding the exponential of the Hamiltonian 
in Chebyshev polynomials, and exploiting the trace op- 
eration performed when evaluating the expectation value 
of an observable, it is possible to reduce the evolution of 
any observable to a one-dimension function that can be 
evaluated directly. 

We have also presented an optimal algorithm to per- 
form such a calculation, and we have shown how this 
new method can easily compete in term of computa- 
tional costs with a variety of model reduction approaches, 
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FIG. 5. t against \f(t)\ plot, the dot is the value at time St. 



whilst maintaining the approximation errors below a cho- 
sen threshold. 

Moreover we have also studied the Zero Track Elimi- 
nation method that has appeared in the literature. We 
have demonstrated that the ZTE method cannot be a 
reliable method for every system, but can be a powerful 
alternative in the restricted setting of NMR simulations. 

G.M. is very grateful to Arieh Iserles for useful sugges- 
tions at the start of this work. 



Appendix A: The DEC algorithm 

We provide here a detailed description of the algo- 
rithm. 

inputs: L hermitian matrix n x n, go vector 

outputs: expectation value f(t) evaluated at f(jAt), 

3 I V. 

1. evaluate a,/3 via Lanczos s.t. a < cr(L) < (3 

2. scale L and get L s 

3. evaluate To = Idgo, T\ = L s qq 

4. while ||cfe|| < e 

5. Cfc = (2 — 5k^){—i) k Jk(rS), t = total time 
6- Tk+\ — 2TfeL s — Tfc_i 

7. store f k = Trace{(T fe £> )<2} 

8. end 

9. for j = 1 : N 



10. re-evajuafce the c k at different time t = jdt 

11. /(j) = £c fe T fe 

12. end 

Appendix B: A counter-proof of the ZTE elimination 

We provide a specific counter-example for the Zero 
Track Elimination. In FigJS] we have plot the absolute 
value of the function: 



/(«) 



-i2nt 



-i2ir(1.001)t 



1 



e -i27r(0.999)t^ (g-g 



If we look at the St evaluated looking at the lowest iso- 
lated frequency we have that 5t ~ 1, within this time 
max\f(t)\ ~ 2 x 10~ 5 but clearly the maximum ampli- 
tude of this periodic function will be 2. 
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